Estimation of Atmosphere and Object Properties in Hyperspectral Longwave Infrared Data

We present a method for atmospheric estimation in hyperspectral longwave infrared (LWIR) data. The method also involves the estimation of object parameters (temperature and emissivity) under the restriction that the emissivity is constant for all wavelengths. The method is analyzed with respect to its sensitivity to noise and number of spectral bands. Simulations with synthetic signatures and signatures from real vegetation are performed to validate the analysis. The proposed method allows estimation with as few as 10–20 spectral bands at moderate noise levels. Using more than 20 bands does not improve the estimates. 1.0 INTRODUCTION Multiand hyperspectral image exploitation is a growing field both in remote sensing within the civilian community and in defence applications such as reconnaissance and surveillance. Multiand hyperspectral electro-optical sensors (cameras) sample the incoming radiation at several (multispectral sensors) or many (hyperspectral sensors) different wavelength bands. Compared to a consumer camera that, typically, uses three wavelength bands, corresponding to the red, green and blue colours, hyperspectral sensors sample the scene at a large number of wavelengths (spectral bands), often several hundred. Moreover, these spectral bands can be beyond the visible range, i.e., in the infrared domain. Each pixel thus forms a vector of measurements from the different bands. This vector, the observed spectral signature, contains information on the material(s) present in the scene, and can be exploited for detection, classification, and recognition. At some wavelengths, the influence of the atmosphere between the observed object and the sensor is significant, and it is often necessary to perform atmospheric correction before further processing. 1 Effects of the atmosphere are typically known to some degree; there are standard simulation tools that provide good approximations. This paper treats methods for refining these approximations from observed spectral signatures of unknown objects, i.e., in-scene atmospheric estimation. 1.1. Scenario A quite specific scenario is considered in this paper, which leads to several assumptions in the treated estimation problems. In our scenario, a reconnaissance aircraft equipped with a nadir-looking sensor has passed over an area, and we want to analyze the observed objects on the ground. The flight altitude is approximately one kilometer. Since we know the time, date, and location, we also have a quite good idea of the atmosphere parameters: humidity and temperature profiles, aerosol and carbon dioxide contents can be read from standard atmosphere models and adjusted using local measurements (e.g., the aircraft can measure the air temperature). Report Documentation Page Form Approved OMB No. 0704-0188 Public reporting burden for the collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information, including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson Davis Highway, Suite 1204, Arlington VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to a penalty for failing to comply with a collection of information if it does not display a currently valid OMB control number. 1. REPORT DATE OCT 2009 2. REPORT TYPE N/A 3. DATES COVERED 4. TITLE AND SUBTITLE Estimation of Atmosphere and Object Properties in Hyperspectral Longwave Infrared Data 5a. CONTRACT NUMBER 5b. GRANT NUMBER 5c. PROGRAM ELEMENT NUMBER 6. AUTHOR(S) 5d. PROJECT NUMBER 5e. TASK NUMBER 5f. WORK UNIT NUMBER 7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) Swedish Defence Research Agency, FOI P.O. Box 1165 SE-58111 Linköping Sweden 8. PERFORMING ORGANIZATION REPORT NUMBER 9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSOR/MONITOR’S ACRONYM(S) 11. SPONSOR/MONITOR’S REPORT NUMBER(S) 12. DISTRIBUTION/AVAILABILITY STATEMENT Approved for public release, distribution unlimited 13. SUPPLEMENTARY NOTES See also ADB381583. RTO-MP-SET-151 Thermal Hyperspectral Imagery (Imagerie hyperspectrale thermique). Meeting Proceedings of Sensors and Electronics Panel (SET) Specialists Meeting held at the Belgian Royal Military Academy, Brussels, Belgium on 26-27 October 2009., The original document contains color images. 14. ABSTRACT We present a method for atmospheric estimation in hyperspectral longwave infrared (LWIR) data. The method also involves the estimation of object parameters (temperature and emissivity) under the restriction that the emissivity is constant for all wavelengths. The method is analyzed with respect to its sensitivity to noise and number of spectral bands. Simulations with synthetic signatures and signatures from real vegetation are performed to validate the analysis. 15. SUBJECT TERMS 16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF ABSTRACT SAR 18. NUMBER OF PAGES 14 19a. NAME OF RESPONSIBLE PERSON a. REPORT unclassified b. ABSTRACT unclassified c. THIS PAGE unclassified Standard Form 298 (Rev. 8-98) Prescribed by ANSI Std Z39-18 Estimation of Atmosphere and Object Properties in Hyperspectral Longwave Infrared Data 12 2 RTO-MP-SET-151 UNCLASSIFIED/UNLIMITED UNCLASSIFIED/UNLIMITED The sensor is a hyperspectral sensor in the longwave infrared (LWIR) domain with around 200 spectral bands in the 7.5–12 μm interval. The final goal of the problem studied here can be formulated as follows: Given a hyperspectral observation of a scene in the longwave infrared domain, estimate the emissivity spectrum and the temperature of the observed objects. This problem is commonly referred to as temperature-emissivity separation (TES), and we will study a special case where we put heavy restrictions on the emissivity. A first step in the TES process is to be able to correct for the atmospheric influence on the observation. The immediate goal of this study is therefore to estimate relevant atmosphere parameters. This is however not completely separable from the TES problem, and, in fact, for certain special cases, the TES problem is solved simultaneously. The secondary goal is to establish sensor requirements, i.e., how many spectral bands are needed at what noise level. 1.2. Related work Despite a large research effort on atmospheric estimation and correction and its impact on detection, there is little work with the same focus as here, i.e., focusing on minimal sensor requirements. Also, the methods in the literature require either high-resolution spectral data, or use spectral features in other wavelengths than the ones we study here. 2, 3 The most well-known methods for atmospheric correction are probably AAC 4 and ISAC. 5 Extensions include ARTEMISS, 6 where a database of atmospheric transmissions is used and a smooth emissivity spectrum is retrieved. The type of parameterized models for the atmosphere that we use here is also used by Fox et al. 7 and Chandra et al. 8 The impact of atmospheric correction on target detection is studied by Yuen and Bishop. 1 1.3. Outline Section 2 explains how we model the observed objects, the atmosphere, and the sensor. In order to adapt the models to observed data, their parameters are estimated using the optimization procedures described in Section 3. These optimization procedures are analyzed with respect to sensitivity in Section 4 and simulations are performed to validate the analysis and to examine the effects of different spectral resolutions in Section 5. Conclusions are drawn in Section 6. 2. MODELLING In this section, we describe how we model the three important entities present in our scenario: the observed object, the atmosphere, and the sensor. The type of models used is non-controversial, even if other modelling methods exist, for example subspace models 8 and databases. 6 2.1. Overview We want to estimate and study the thermal radiance from an object (or from several objects) represented by a vector Lo or a continuous function Lo (λ). In a laboratory setting, where the sensor is close to the object, the problem would be considerably less complicated. In our scenario, the sensor is airborne around one kilometer above the observed objects, and we thus need to take atmosphere effects into account. We model the atmosphere by a function A(·), and the at-sensor radiance Ls is thus a function Ls (λ) = A(Lo (λ)). Estimation of Atmosphere and Object Properties in Hyperspectral Longwave Infrared Data RTO-MP-SET-151 12 3 UNCLASSIFIED/UNLIMITED UNCLASSIFIED/UNLIMITED Disregarding any spatial effects, the sensor is modelled as a global scaling factor and as a spectral filter for each spectral band. The spectral filtering is the function that samples a continuous function L(λ) and outputs a measurement vector z. The sensor is modelled by a function S(·), and our observation vector is modelled by z = S( Ls(λ) ). Below, some details about sensor, atmosphere, and object models are given. 2.2. Discretization In practice, we perform all computations on wavelength-discrete vectors. The resolution of these vectors are either determined by the sensor’s spectral resolution or by the high-resolution simulations performed by MODTRAN. These simulations give high-resolution transmission/radiance vectors (every inverse centimeter) that for our purposes serve well as (approximately) continuous functions. The blackbody function determining the object radiance can naturally be computed in any resolution. 2.3. Sensor model Recalling that we do not make any spatial considerations here, the sensor can be modelled purely as a spectral function S : L→R. The sensor inputs an observed radiation L(λ) (a continuous function of λ) and outputs a discrete vector y of N measurements. Each measurement yn is a sample of L(λ), and it is sampled using a response function Rn(λ) centered around λn. The response function Rn(λ) is typically modelled as a Gaussian function with a certain width δn. The sampling is thus expressed by a function P : L 2 →R N such that P( L(λ) ) = [P1( L(λ) ), . . . , PN( L(λ) )] T . Unfortunately, the sensor is not perfectly linear, but we will in the following assume that the sensor is calibrated as to be linear and with a global scaling factor s compensating for variation in gain and flight altitude. Thus, we model the measurements as yn = Sn( L(λ) ) = s · Pn( L(λ) ). The discrete sensor function samples high-resolution vectors using the Gaussian spectral filter. Our discrete sensor model, is S(L) = [S1(L) ... SN(L)] T which we can write in matrix form so that S(L) = s P L. 2.4. Atmosphere model The atmosphere between the observed object and the sensor influences the observed signature due to transmission, radiance, and spatial distortion (the latter due to turbulence or scattering into the path). Again, disregarding spatial effects, we model the atmosphere as a spectral function A : L→L. The atmospheric transmission is modelled as a multiplicative transmission in each wavelength, and since the atmosphere is also a source of radiation in itself, the model of the observation is A( L(λ) ) = τ(λ) · L(λ) + L ↑ (λ), where L ↑ is the up-welling path radiance. Note that this model of the atmosphere assumes a homogeneous layer of air between the sensor and the observed object. The atmosphere is variable in temperature and humidity (the water vapour content), and these parameters influence the atmospheric transmission and radiance. Assuming that the atmospheric transmission τd(λ) Estimation of Atmosphere and Object Properties in Hyperspectral Longwave Infrared Data 12 4 RTO-MP-SET-151 UNCLASSIFIED/UNLIMITED UNCLASSIFIED/UNLIMITED without water vapour and the transmission of (a certain amount of ) water vapour τv(λ) are known, we can parameterize the atmospheric transmission as τ(w,λ) = τd(λ) · τv(λ) w , (2) where the humidity scaling w is the relative amount of water vapour, i.e., w = v/v0, where v is the actual amount of water vapour and v0 is the amount of water vapour giving the transmission τv(λ). The atmosphere radiance depends on the atmospheric transmission and the temperature profile ta(z), i.e., the air temperature at each altitude z. Here, we simplify the temperature profile to a constant ta and thus use a one-layer atmosphere model. The emissivity spectrum of the atmosphere is c(λ) = 1 − τ(λ), and the radiance is thus expressed as L ↑ (w, ta, λ) = (1 − τ (w,λ)) · Lbb (ta, λ) (3) In the discrete case we have τ(x) = T(w) x, where T(w) = diag( τ(w) ) is a diagonal matrix with elements 0 ≤ τi ≤ 1. τ(w, λ) depends on the parameter w, as described above. The discrete atmosphere model is thus A(L) = T(w) L + L ↑ = T(w) L + (I − T(w)) Lbb(ta), where ta is the temperature of the atmosphere and Lbb(t) is the radiation from a blackbody with temperature t. Note that the equality (2) becomes an approximation in the discrete case. As long as w is close to one, the approximation is close to equality. The atmospheric transmission used in the experiments in this paper is a simulation of the transmission from ground level to one kilometer above, at midday of a Swedish summer’s day. From this, the upwelling radiance can be computed according to (3). The down-welling radiance L ↓ (λ) is simulated for the same time. 2.5. Object model The radiation at the wavelength λ from an object is the sum of the emitted and the reflected radiation. We write this as Lo(λ) = Ls(λ) + Lr(λ). The emitted radiation is characterized by the object’s emissivity spectrum c(λ) and its temperature to as Ls(λ) = c(λ) · Lbb(to , λ). In the following, we will consider the special case where the observed object is a graybody, i.e., c(λ) = c is constant for all wavelengths. The reasons for using this model are simplicity (two parameters only) and that we in a natural environment have a large chance of finding several graybody pixels. The reflected radiation depends on the incoming radiation E(λ) and the emissivity as Lr(λ) = (1 − c) E(λ) / π. Note that we neglect any directional effects here, i.e., we assume that the incoming radiation and/or the object surface is 100% diffuse. Returning to our scenario with an airborne, nadir-looking sensor, and also assuming that we have pure pixels, the discrete object radiance model at pixel coordinate u = (u,v) is Lo(u) = c(u) Lbb(to(u)) + (1− c(u))L ↓ . where L ↓ is the down-welling atmosphere radiance. Estimation of Atmosphere and Object Properties in Hyperspectral Longwave Infrared Data RTO-MP-SET-151 12 5 UNCLASSIFIED/UNLIMITED UNCLASSIFIED/UNLIMITED x 2.6. Putting it together Summarizing the models above we can model the observed vector z at image coordinate u as z = S( A( Lo(u) ) ) = s P A( Lo(u) ) = s P ( T(w) Lo(u) + (I – T(w) Lbb(ta) ), where T(w) = diag( τd ) diag( τv w ) and Lo(u) = c(u) Lbb(to(u)) + (1 – c(u)) L ↓ • τd, τv, and L ↓ are computed using MODTRAN. • The sensor is characterized by the parameters s, {δn}, and {λn}. • The atmosphere is parameterized by w and ta. • The object radiance is characterized by c(u), to(u), and L ↓ . Given the parameters x = {s, ta, w, to, c} we have a generative model for the observation. Thus, we can write y as a function of a parameter vector x, i.e., y(x). 2.7. Noise model In practice, the measured vector will not equal y exactly, since we have noise and model errors. We have model errors since none of the assumptions above are 100% true, and we have noise added by the sensor. We will model the noise as a being zero-mean, white, additive and with a probability distribution function monotonically decreasing with its square, for example the Gaussian distribution. Thus, our measurement will be z = y(x) + n, where n ∼ N(0,σ). 3. OPTIMIZATION PROCEDURE Minimizing the residual between an observation and a parameterized object-atmosphere-sensor model is an optimization problem in several variables. The problem is non-linear in some of the parameters and linear in some. Having a model vector y that depends on a parameter vector x and a measurement vector z, we estimate the parameters by minimizing f(z,x) = ||k(z,x)|| 2 = ||z – y(x)|| 2 with respect to x. The objective function f is non-linear and a suitable optimization algorithm should be applied to find the optimal x. The problem is simplified by observing that c and s are linear parameters that can be solved for in a linear sub-problem given values on to, ta and w. Depending on our assumption about the object, this is done in slightly different ways. Using our previous assumptions on the noise, we solve in a least-squared sense, and our notation for a constrained linear least-squares problem is to find